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ABSTRACT 


The viscous characteristic analysis for supersonic chemically reacting flows 
as reported in NASA CR 111783, has been extended to include provisions for 
analyzing embedded subsonic regions. This report describes the numerical 
method developed to analyze this mixed subsonic-supersonic flow fields. A 
discussion of the boundary conditions related to the supersonic-subsonic 
and subsonic-supersonic transition as well as a heuristic description of 
several other numerical schemes for analyzing this problem is included. An 
analysis of shock waves generated either by pressure mismatch between the 
injected fluid and surrounding flow or by chemical heat release is described 
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I. INTRODUCTION 

One of the interesting phenomena produced by combustion in a supersonic flow . 
is the possibility of producing transition from supersonic to subsonic flow 
due to the heat addition. This transition has been analyzed in the past for 
the case of one dimensional flow, however, a more complex flow field is usu- 
ally generated in a supersonic combustion process controlled by mixing. In 
the flame zone of this flow, where the static temperature of the gas in- 
creases, the local Mach number can decrease, even if the static pressure de- 
creases, because the value of the local speed of sound increases. Therefore, 
a transition from supersonic to subsonic flow can take place in a localized 
region of the flow where the pressure changes slowly and the region can be 
completely surrounded by supersonic flow. 

The transition from subsonic to supersonic flow can occur through several 
mechanisms; when the flow surrounding the subsonic region is cooler than the 
\ flow in this region (because combustion takes place only in a limited region 
of the flow), then, because of mixing of the combusted gases with the surround- 
. ing air, the temperature of the gas in the 'subsonic region gradually decreases 
and the flow can become supersonic without pressure variation, or even with a 
small static pressure raise. The second possibility is that the stagnation 
temperature remains roughly constant but the pressure decreases, the flow 
accelerates and again becomes supersonic. The subsonic flow field region is 
essentially the opposite of the classical transonic region about airfoils, 
where the flow becomes locally supersonic, because the pressure decreases and 
then again becomes subsonic downstream, usually through a shock, and a local- 
ized region of supersonic flow is embedded in a subsonic stream. Extensive 
discussions on the existence of transonic flows without shocks have taken 
place in the past for this type of flow. It is now well recognized that smooth 
solutions (without shocks) are special solutions of this transonic problem and 
require special boundary conditions. 

The possibility of transition from supersonic to subsonic flow due to com- 
bustion has been recognized in the past. In Reference (1) this possibility 
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was discussed qualitatively and some of the important aspects of the presence 
of such subsonic regions were described, however, a detailed treatment of 
this problem is still lacking. A quantitative understanding of the physical 
aspects of such a transition are of primary importance for the development of 
a practical scramjet engine, in view of the fact that the efficiency of the 
engine, at flight Mach numbers of the order of 6 to 8, depends strongly on 
the value of the static pressure and local Mach number at which the heat ad- 
dition takes place. High efficiency and high specific impulse requires a 
low local Mach number where heat is released due to chemical reaction. If 
the Mach number before combustion is low, during the combustion process sub- 
sonic regions are formed, even if the pressure locally does not change be- 
cause of the increase of the local speed of sound. In addition, the scram- 
jet must be able to operate at lower flight Mach number than design. For 
these conditions, large regions of the burner have subsonic flow. 

ATL has developed for NASA in the past, methods where the flow fields of or 
supersonic streams with chemical reactions can be analyzed (Reference 2). 

Such numerical methods permit analyzing the flow between discontinuities, 
provided that the entire flow remains supersonic, in view of the fact that 
such a method uses "viscous characteristics" as first described in References 
(3) and (4). In the present phase of the work, a numerical method capable of 
determining the formation of combustion shocks has been developed. In addition, 
the possibility of numerically determining the region of transition from 
supersonic to subsonic flow, and then subsonic to supersonic as can occur 
in the flame has been investigated. 

The numerical investigation of a transonic problem, as described requires the 
development of a numerical scheme, valid in the transonic region, that can be 
coupled with the "viscous characteristic program". In addition, it requires 
a clear understanding of appropriate boundary conditions for the subsonic 
regions. A problem similar to the problem related to the existence of smooth 
subsonic-supersonic transition over a two dimensional profile near M=l, exists 
for this flow; therefore, the first step in this effort is to generate a clear 
understanding of the boundary conditions related to the transition problem. 
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II. BOUNDARY CONDITIONS RELATED TO THE TRANSITION FROM SUBSONIC TO ' ' 

SUPERSONIC FLOW AND FROM SUPERSONIC TO SUBSONIC FLOU 

In the combustor of a scramjet engine, the flow downstream of the burner is 
supersonic, and the flow reaching the burner is also supersonic, therefore, 
any subsonic region generated locally by combustion is contained between 
these two supersonic regions. Consider first the case in which the subsonic 
region is completely imbedded in a supersonic stream, as shown in Figure (1), 



The hydrogen from the injector mixes with the air outside and combustion 
takes place; the pressure rises because of the decrease in density of the 
stream due to combustion; the speed of sound increases, and in this high 
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temperature region the flow can become subsonic. Since the duct in which 
burning occurs diverges and the downstream pressure is much lower than the 
pressure at the injector face, the flow must again cross a sonic line. 

If the temperature in the high temperature zone is higher than the average 
temperature after combustion (case of 4><1), the flow outside the high tem- 
perature region will remain supersonic. Therefore, the shape of the sonic 
line is as shown in Figures(l) and (2). The static pressure distribution 
between A and B does not uniquely define the flow and the line M=1 is not 
a line of constant static pressure temperature; and velocity as in the case 
of the wing, because the stagnation temperature changes from point to 
point. Hence, several different pressure distributions are possible de- 
pending on the amount of diffusion and chemical reaction taking place be- 
tween A and B. The pressure can first increase and then decrease, can 
remain approximately constant, can continuously decrease, or can slowly in- 
crease. In this last case, the variation of the 6 speed of sound due to dif- 
fusion must be larger than in the other cases. 


SONIC LINE 



FIGURE 2. 
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It is not clear & priori, if the presence of a closed sonic line in the 
flow, as shown in Figure (2) is physically possible. The boundary conditions 
required in. order to obtain such smooth conditions must be investigated first. 

In order .to understand the physics of the problem it is convenient to analyze 
the problem in two steps. The first step is related to the transition from 
subsonic to supersonic velocity. Then we can transform the problem in such 
a way that the properties of the subsonic region can be controlled indepen- 
dently of the supersonic region, as shown in Figure (3). We assume that 
flow 1 can be controlled independently of flow 2 and that flow 1 is ini- 
tially subsonic. 



FIGURE 3. 
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The simplest case of such a flow occurs when stream 1 is uniform and initially 
subsonic (one dimensional) at the exit of the jet and stream 2 is supersonic. 
If we neglect transport properties and assume both flows to be inviscid, the 
flow properties as specified along the line, AC yield a unique relation be- 
tween the pressure and the angle of the streamline AB, where CB is a charac- 
teristic line emenating from C. The unique relation between p and 6 cor- 
responds to a unique relation between p(x) and y(x) since 

x 

y = y tan 8 dx 
o 

A similar relation can be determined for streamline AB. For example, if 
we assume that the pressure in region 1 is independent of y (one dimensional 
flow), we have a simple relation between y and p given by the continuity 
equation. Thus, a step by step calculation can be performed where at each 
step Ax between A and B, the variation Ap in the step can be assumed as a 
parameter and a single calculation can be performed (Figure 4). The equa- 



FIGURE 4. 
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tion of characteristics applied along Cb yields a relation between p b and 


For a two dimensional flow 


where 


2 (Pb-Pc) _ (6b -6c) = ' Q 

Y (Pb+P c ) . cosu . siny 


tanp= j (tany b + tanu c ) 


(i) 


Then, for a given value of Ap (i.e., p b ~p a ), P5 is given and e b is calcula- 
ted. Then, Ay is determined from 


2 (tan0 a + tan9 b ) Ax 


The one dimensional relation yields 


yb. = (El) 

ya Pb 


Pa, ^ 


III 
P a Y, 

1 - (jr> 

1 H 0 


111 

\i-(V 


1/2 


1 / 


where p Q is the stagnation pressure of stream 1 


( 2 ) 


(3) 


Only one value of A p satisfies Equations (1), (2) and (3). Then^the varia- 
tion of y as a function of x can be obtained. However, in the region where 

M, becomes equal to one ^ =0. Then. a solution can not be found un- 

Qp M=1 

less 0 0 and tends to zero for the external flow in the region where 

+ 1. If this does not occur, we have a "choking' 1 condition and the ini- 
tial subsonic flow distribution assumed in Figure (3) is not physical. 
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A "choking' 1 condition can be found independently of the approximation used 
in the analysis of the subsonic flow. The physics of the simple case shown 
in Figure (3) is clear. The pressure at station A has been assumed arbi- 
trarily. If choking occurs, this implies that an incorrect value of the 
pressure at A has been assumed. If we change the value of at A and, 
therefore, the corresponding value of p^, the shape of the streamline AB 
changes. If the value of 0 of the streamline AB at the choking point is 
negative, then the mass flow of stream 1 must be decreased; hence p A 
a given p Q must increase. This change tends to increase the value of 9 at 
the choking region. Therefore, by an iterative process the value of p^ 
that gives a physical solution can be determined. The problem is sub- 
stantially more complex if flow 1 is not assumed to be one dimensional, 
and transport properties and chemistry are included in the analysis, as 
will be discussed later on. However, the controlling mechanism is the same. 

For any given initial flow distribution in the supersonic region, and given 
stagnation conditions and channel shape of the subsonic flow, a single solu- 
tion can always be found that permits the subsonic flow to cross the M=1 line, 
which corresponds to a given value of the static pressure at A. Therefore, 
the correct boundary conditions for the problem requires a selection of the 
pressure at A that avoids "choking"^ and the assumption that the downstream 
pressure is sufficiently low. Then ; the initial flow properties along AC 
and the geometry defines the problem. 

Let us now consider the case depicted in Figures (1) and (2). In this case 
the flow is initially supersonic, therefore, the flow field in front of the 
sonic line is completely determined by the initial conditions. In Figure 
(5) the flow along the line a is given. Then, the flow along line b is 
uniquely determined. If the flow becomes subsonic due to chemical reaction, 
the flow properties along the sonic line S are completely determined; therefore, 
only one solution can be found that gives a smooth transition from supersonic to 
subsonic flow and a unique flow can be determined along line S where M=I hav- 
ing an embedded smooth subsonic region. Hence, the possibility of chang- 
ing a parameter, equivalent to changing the value of the pressure p^ of Figure 




(3), (in the case that choking occurs because of the downstream conditions) 
is not available for this flow. However, if downstream, the flow Mach num- 
ber tends to increase and the subsonic region tends to decrease, then the 
sonic line crosses all the streamlines and closes as shown in Figure (6), 



FIGURE 6. 
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Therefore, the possibility of having the choking condition described be- • 
fore, still "exists. For some conditions, the flow cannot undergo transition 
from subsonic to supersonic because the flow is choked; therefore, some modi- 
fications must be introduced upstream to change the initial conditons of the 
subsonic flow. The only existing mechanism for this adjustment is the forma- 
tion of a shock. The choking condition produces a disturbance that moves up- 
stream of the sonic line (Figure 7) forming a shock. The transition from 
supersonic to subsonic then occurs by means of this discontinuity as shown in 
Figure (7). 

The position and shape of the shock depends on the choking conditions, which 
are dependent on the flow properties along AB. Then, in the general case, 
for any given flow distribution along AB, the transition from supersonic to 
subsonic occurs through a strong shock (subsonic flow behind), the position 
of this shock being highly dependent on the flow properties in the super- 
sonic region outside the choking region (region CC of Figure 7). 

Similarly to the transonic case, a smooth transition can exist for specially 
selected boundary conditions and can be calculated by specifying only a part 
of the flow properties along the initial station AB, for example, along AD 
of Figure (7). Then, the flow properties along DB can be determined by an 
inverse process where the transonic region CC is determined first and the 
region DCCB is determined later on. 

The determination of the shock cannot be obtained by a direct (marching) 
calculation as it depends on the downstream flow conditions; it should be 
obtained by a procedure similar to that used for transonic flow analysis. 

The procedure that appears the most feasible for this problem is an itera- 
tion procedure where a smooth transition from supersonic to subsonic is 
assumed first and a ‘'choking" error is determined at the subsonic to super- 
sonic transition in terms of the amount of mass flow that cannot cross the 
region of M=l. Then, a shock position, s, is assumed ahead of the first 
transition. This position, s, defines a new flow field . Such a flow 
field defines a new error associated with choking at the second transition. 



M >1 


SONIC LINE 



FIGURE 7. 
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The error obtained is smaller than the previous error. A program that can 
perform this type of analysis and perform these iterations can be developed, 
however, it is not part of this effort. The effort supported by NASA in the 
mentioned contract is limited to the development of a program that analyzes 
only the inverse problem. This analysis is described in detail in the follow 
ing sections. 
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III. GOVERNING EQUATIONS 

The analysis developed for calculating the subsonic region of the flow must 
be consistent with the "viscous characteristic" method employed in the 
supersonic region and must be capable of analyzing a viscous, reacting flow 
with both lateral and transverse pressure gradients. The pressure p and 
flow inclination e have been selected as independent variables since the 
equations written in terms of p and 0 do not explicitly contain the entropy 
as a flow variable, which can have extremely large gradients in the combus- 
tion region. Use of equations written in terms of velocity components, 
could lead to numerical difficulties regarding mesh size in the combustion 
zone since the entropy gradients would have to be evaluated by means of 
finite difference formulas. 

The governing equations are the well known "viscous-inviscid" equations em- 
ployed in higher order boundary layer and viscous flow field analyses 
(References 5 and 6) with the finite rate chemistry terms included. The 
equations are as follows: 

Global Continuity : 

^ +pqf^+^ a sTn6 =0 14) 

(J=0 for two dimensional fow and J=1 for axi symmetric flow) 


S-Momentum: 


pq 


M+ M. = 


3S 3S 


= s? ^ + y cos 6 11 frl = 


N-Momentum: 


2 36 ap _ - 

pq 3? + an 0 


(5) 


( 6 ) 
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Energy : 


C p q — - q = ? 

"p pq as Re ( Y -i)M‘ 


00 00 


^ 3T + |,L e 3T _ r 3a i 

r £r- T-r + TTT- rr ^ t. 


3 ^ II + Iras B 

¥?T PF sn + 7 cos 6 PF aW ' Pr an “ “p. an 




2-, 


- E *i h i = S 2 ' Z *i h i 


(7) 


Species Conservation : 


pq 


3a • i 

1 = J_ 
3s Re 


an Pr an y Pr cos 6 w an 


( 8 ) 


+ to. = S- + to. 
1 m 1 


State: 


P = 


W pT 

OCT 

y m 2 w 

' OQ 05 


where 


a i 

W = s — 

m. 


-1 


0) 


These equations are written in an intrinsic coordinate system with s along 
and n normal to the streamlines. We have assumed that transport effects 
are produced only by gradients normal to streamlines. 
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IV. HEURISTIC DESCRIPTION OF SEVERAL NUMERICAL SCHEMES FOR ANALYZING, 

SUBSONIC ZONES 

A marching technique (as described in Section II) may be developed using a 
one dimensional approximation for the subsonic stream. Such an analysis 
is described in Reference (7) for analyzing supersonic air-ejector flow 
fields. In using a two dimensional description of the subsonic stream,, 
the mathematical nature of the problem changes. In this representation, 
the normal momentum is introduced, an equation not accounted for in the 
one dimensional approximation. The governing flow equations in the super- 
sonic region of the flow field are hyperbolic-parabolic in nature and hence 
may be solved by a marching technique; but, in the subsonic portion of the 
flow field, the flow equations are elliptic. This situation presents 
significant difficulties in attempting to analyze viscous, combusting flow 
fields, since .both mixing and chemistry are analyzed numerically by marching 
along the flow streamlines. 

While a marching scheme for an elliptic system yields an improperly posed 
problem mathematically, it has provided solutions of the inverse blunt body 
problem. Hence, a marching type numerical method was envisioned as a pos- 
sible approach to the solution of this problem. Regarding such a numerical 
approach, (Reference l ) , “fundamental questions arise with respect to the 
uniqueness and existence of a solution and with respect to stability and con- 
vergence of calculated procedures." Then, the possibility of obtaining 
physical solutions to elliptic problems by a marching scheme is highly de- 
pendent on the scheme utilized. 

A marching scheme was developed employing the governing equations in finite 
difference and using the numerical scheme described in Appendix I. This 
scheme for analyzing the subsonic region was incorporated into the "viscous- 
characteristic" program and test cases were performed to analyze mixed subsonic 
supersonic flow fields. Both "direct" problems (where an upper subsonic 
boundary shape was iterated upon to yield a lower subsonic boundary that 
satisfies either a wall, axis or characteristic compatibility constraint) 
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and "inverse" problems (where the upper subsonic boundary shape is pre- 
scribed, yielding a wall shape at the lower subsonic boundary) were at- 
tempted. In all cases that were run (both of the "direct" and "inverse" 
type) oscillations developed in both the pressure and flow deflection 
profiles after marching several axial stations which grew in magnitude 
and caused the program to terminate. It was felt that a polynomial repre- 
sentations for the pressure and flow deflection in the subsonic region 
might alleviate the instabilities obtained with the finite difference ap- 
proach. This polynomial scheme is described in Appendix II. This 
method permits the calculation of a mixed subsonic-supersonic flow fields 
by a marching technique employing "viscous-characteristics" in the super- 
sonic portion of the flow field and a multi-strip integral technique for 
the subsonic portion. A simple two dimensional test case was run employ- 
ing the initial profiles shown in Figure (8), with a wall producing a 
rapid expansion imposed as an upper boundary, as sketched in this figure. 

A solution was sought that would accelerate the flow producing an all 
supersonic flow field downstream of the initial station. 

For this test case, the medium was air, and the viscosity was set to zero, 
to simplify understanding the elemental physics involved. The flow is, 
however, rotational and non-homentropic. The matching point selected on 
the initial profile was at a Mach number of 1.03. 

The marching scheme, as described in Appendix. II, entails an iteration for 
the pressure at the matching point, the correct pressure being the one that 
passes the appropriate value of mass flow. Cases 1, 2 and 3 are distinguished 
only by the differences in the initial flow deflection (e) profiles. As 
indicated by streamwise pressure and Mach number variations at the axis 
Figure (9), the flow solutions differed substantially for these cases. In 
case (1), the flow accelerated to a local minimum section which was ap- 
parently larger than the critical area required to accelerate the flow to 
the supersonic branch, hence the flow decelerated downstream of this station. 
In case (2), the flow smoothly accelerated from subsonic to supersonic, while 
in case (3) the flow reached a station where local choking occurred. Hence, 
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only case (2) could be analyzed by a direct marching approach. Case (1) . 
cannot be handled by a direct marching technique since the downstream 
conditions are not known a priori while case (3) requires the addition 
of shock waves to adjust the incoming mass flow. Figure (lp) indicates 
M, P and e profiles at several axial stations for case (2). 

When the polynomial method was applied to the analysis of viscous-combusting 
flow fields the following difficulties were encountered: 

(1) Both the pressure and flow deflections exhibited oscillations 
at the upper and lower bounding subsonic streamlines which 
tended to rapidly grow in magnitude. It was found that 
these oscillations were strongly a function of the size of the 
marching step in the streamwise direction; the larger the 
axial step size, the smaller the amplitude of the oscilla- 
tion. This is analogous to the numerical difficulties en- 
countered in the blunt body problem, but cannot be simply 
resolved since both the chemistry and the mixing impose 
limitations on the allowing axial step size. 

(2) When one of the bounding streamlines was a wall, it was 
found that the flow deflection polynomial e (y ) could not 
physically describe a subsonic-supersonic or supersonic to 
subsonic transition at this wa 11. 

(3) In attempting to run problems of the "direct" type, local 
choking consistently occurred causing the program to termi- 
nate. 

On the basis of the results obtained with both numerical schemes developed 
the following conclusions were reached: 

(a) A marching scheme employing a finite difference approxi- 
mation of the governing equations (Equations 4 - 9) of 
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all grid points appears to be unstable for both * 
direct and indirect problems. 

(b) A marching scheme employing polynomial expressions 
for the pressure and flow deflection yields oscilla- 
tions at the bounding streamlines; cannot represent 
bounded flow fields with subsonic flow adjacent to 
one of the boundaries and yields local choking for 
direct problems. 

The numerical scheme described in the following section yields a workable 
method for problems of the inverse type, for embedded subsonic regions ad- 
jacent to walls. 
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V. NUMERICAL PROCEDURE FOR THE ANALYSIS OF EMBEDDED SUBSONIC 

REGIONS 

A numerical method of the inverse type has been developed for the analysis 
of embedded subsonic regions adjacent to walls. A typical flow field is 
depicted in Figure (11). 



/ 
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The flow field contained between Station (1) (initial profile) and Station 
(2) is entirely supersonic and hence may be analyzed by “viscous- 
characteristics" as described in Reference (2). It should be pointed out 
that the "viscous-characteristic" program as described in Reference (2) 
was catered to analyzing expanding supersonic nozzle flow fields. In this 
current effort, combustion generated compression fields are being analyzed, 
hence major changes have been incorporated into the basic "viscous- 
characteristic program. These changes include the capability of calculat- 
ing the initial mixing region including the determination of an initial 
underexpansion shock; carrying shocks as discrete discontinuities in the 
flow field; the detection and calculation of envelope shocks produced by 
the coalescence of compression waves; and grid spacing logic for analyzing 
a mixing region of small transverse extent in a nonuniform supersonic flow 
field. These revisions will be described in later sections of this report. 


Consider the flow field at Station (2), where the program has first en- 
countered a grid point with a Mach number M <_ 1,01. At this station, the 
lower wall (AB) becomes the lower subsonic boundary and a slightly super- 
sonic streamline (F6 with M ^ 1.1 to 1.2) becomes the upper subsonic bound- 
ary, as depicted in Figure (12) 



Sonic Line 


Lower Subsonic 
Boundary 


FIGURE 12. 
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The shape of the upper subsonic bounding streamline (F-G) is specified as • 
the polynomial 

“2 —3 

Y fg = A + B X + C + D g~ + (10) 

where X = X - X p 

The first four coefficients in this polynomial are related to the streamline 
shape at F as follows: 


B = (jg) = tane F 
r 

.2 6. 

c -0 = (— §-> 
dx p cos e p 

ri 3 e + 3tane (e ) 2 

D - (H) - — 4 1 

dx p cos e 


(lla) 

(lib) 

(He) 

(lid) 


The quantities e s and e can be obtained from the y derivatives of flow 
quantities at F. Consider the following system of equations; 


/ (M 2 -l) 


0 

sine 



0 0 (yPM 2 )^ 


1 p s! 


,A \ 

1 (yPM 2 ) 0 


P n 


0 

cose 0 0 j 


0 s 


I p y 

0 sine cose / 

0 1 


1 Q 

/ 

n 1 


\ y I 


(12) 


which are modified continuity equation, normal momentum equation and the 
geometric relation a/ay = cose a/an + sine a/as applied to both P and e. 
The righthand side terms may be evaluated at F since they consist of dif- 
ferentiation in the y direction only. (A is related to the mixing and 
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chemistry terms as described by Equation (AI-3). A 'derivation of the 
modified continuity equation may be found in Reference 2). Solution of 
this system of equations then yields e $ , P s> e n and P^ at point F. Ex- 
pressions for 0 s and P g are given by Equations (AI-1) and (AI-2) in 
Appendix I. Then, . 


e 


n 


A - (M 2 -l) P s 
yPM 2 


(13) 


and 

p = . yPmV 
n s 


(14) 


Consider 3/8n of the modified continuity equation, a/as of the normal mo- 
mentum equation and a/ay by both e g and e This yields the following 
system of equations: 



0 0 

0 (yPM 2 ) 

cose sine 
sine 0 



where 

a = A n - 2MM n P s - (yPM 2 > n 6 n 
b = - ( Y PM 2 ) S 8 S 
C = - (ej 

y 



This system may be solved for e ss yielding 


(15) 


(16a) 

(16b) 

(16c) 

(16d) 
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o 

6 ss = (^ n cos8+ 2 cos9M ^ n P s + (YP M ) n COS 0 { 0 n ) 

(-yPM 2 ) cos 2 e(M 2 -l)e +( Y PM 2 )sine(e ) (17) 

s y 

{-yPM 2 COS0 (0 n ) )/ (yPM 2 [M 2 COS 2 0-1] ) 

Note that a/ay of e and e are obtained numerically by obtaining these 

w 0 I 

quantities at grid points above and below F, while A n is approximated by 
consistent with the assumption in the governing equations that for the 
transport terms 3/an ^ a/ay. 

With the upper boundary streamline now specified (if we do not impose higher 
order terms), the pressure distribution could be obtained along this stream- 
line independent of the subsonic flow field, if there were no viscous effects 
Due to the presence of viscosity, the following iterative procedure is used 
to calculate the combined flow fields. 


The lower wall is prescribed by the polynomial 



= a + b x + 




( 18 ) 


where a, b and c are prescribed by the relations (11a), (11b) and (11c) pre- 
scribed at point A. The term d is the parameter iterated upon in the problem 

if a value of d is assumed, the lower wall AB is specified. The pressure 

distribution in the y direction is prescribed by a parabola 


P (y) - a + Sy + cy (19) 

and initially (when first encountering a subsonic zone), the existing pres- 
sure distribution must be fit with this parabola to avoid spurious pressure 
gradients in the subsonic region due to "wiggles" in the actual pressure 
distribution. 


In determining the flow properties at the next axial station, the axial step 
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ax taken is governed by the allowable stability criterion described in 
Appendix III), the supersonic portion of the flow field (points above F ' ) 
may be calculated by "viscous-characteristics". With the flow deflection 
known at point F', the pressure is obtained using the compatibility relation 


B 1 ^ P F'“ P K^ “ ( 0 F‘~V + ^ S y n0 p. K + B 2 + B 3 ] B 4 ax=0 
applied along F'K where 


d — / si nycosy \ 

B 1 ( yP > 


_ Sj (y-l)S 


B 2 2 yPq 
pq 1 M 


2 W 6 \ 

pq i m i 


B, = - ti- (— ) + W i (li- X) j 

3 T 3s chem i 35 ra i . 

chem 




sinp 


4 ” cos(e-y) 


(20) 

( 21 ) 

( 22 ) 

(23) 

(24) 


Note that the chemistry is uncoupled as explained in Reference (2) and the 
term B 3 may be evaluated along the flow streamlines prior to the computations 
of the fluid mechanical properties. 


With the pressure gradient and streamline location known, the velocity, tem- 
perature and species mass fraction may be obtained at F 1 (or any mesh point I) 
by using an explicit finite difference formulation of Equations (5), (7) and 
( 8 ). 


q I q I pjq^pjqj 


2(P I - P I } 2 S 1 AS 




(25) 


where is approximated by derivatives in the y direction 

2 

C - rUL 4. 1 COS0 19.4 /Da 

Sj -u T t~2 + J — 3y I/Re 

ay 


(26) 
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and the derivatives and ^-3- are evaluated at the grid point I by the 

ay 3y 2 

centered finite difference formulas 


non- 


(19.) 

' 8y ' 


n (^1 Ay 2 \ „ 

Uw Alf / ^ 


Ay, 


1+1 Ay 2 H 1 v Ay 2 Ay M I-l. &y x 
_____ 


(27) 


and 


Ay. ay 2 

(A) = 2 [q I»l gr + Ay 2 ' q I + q I-l &y 1 Ay,, 1 

ay 2 j 


^2 


(28) 


where ^ = y t - y^ 
and iy 2 = y I+J - 

The temperature T is determined by the finite difference relation 


T i - T i + As + 


chem 


2 (Pj - P.) 2 S 2 as 

C Pl P I + C Pj P I + ( P qc p ) + ( p qC p )j 


where 


(29) 


S l s (II) 2 + J cose <£(31) + hi ( 3 I )ec !!l 

2 U T L Pr 7T FF FT K zy } 0 y FF ^ay' FF v 3.y ; S, 3y 
ay 


ay' Pi ay 

2 (30) 

+ (y.-D'Mf (|3-) 1/ [Re ( Y _-l) M?] 


CO oo 


All first and second derivatives are evaluated as described by Equations (27) 
and (28). 


_.th 


The i species mass fraction a . is determined by the finite difference rela- 
tion 

/ 8c *i\ 2 S 3 1 As . 

a ij a ij ^ 3S ^ j AS + (p.Tqj+pTqy) 

chem 11 


(31) 
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where 




32a 



+ J 


cose 

y 


L 3a. 
Pr 


(32) 


The remaining properties are made as follows. The average molecular weight 
is obtained from 


a. -1 

W = U i^-)] 
i 

hence the mixtures gas constant is R 
density is 

W„ y M 2 

P 1 00 CO 

p = T “W 


0 

W 1 


(33) 


From the equation of state, the 


(34) 


The specific heat (frozen) is expressed by 

aMT), 

C p “ E C Pi (T ) “i " E (a i lf- } (35) 

where C (T) is a specified function of temperature, thermodynamic data 
^i 

being tabulated from Reference (15). .. 


The ratio of specific heats is 
- C P 

yf ■ £ p- r/c p. 

and the Mach number is given by 


M q ; v R 2 

M* = — l-Vl ' 


/f 


yR 


(36) 


(37) 


The derivative P at F* is computed from P and P . where P is evaluated 
y n s n 

using the normal momentum equation (Equation 14). e s at F 1 is obtained from 
the relation 



3 

cos e, 


[C + D (x F , -Xp )] 


(38) 
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while P s is evaluated by a backward difference. 

The Pjyj polynomial at the new axial station A 1 F * is then obtained employ- 
ing a simple iterative procedure. A value for P at A 1 is assumed (for a 
first guess, the value at A is used), and the polynomial coefficients a, b 
and c are evaluated using Pp , , Py^, and Py^,. This yields the pressure at 
A', hence all flow properties can be evaluated using Equations (25) - (37). 
Then, the normal momentum equation is used to obtain P n and combining this 
with P s an alternate value of P^ is obtained at A‘. The value of P^ ini- 
tially assumed is then perturbed and the process repeated until the two values 
of Py agree to within a specified tolerance. Then, employing the poly- 
nomials at both stations, the pressure gradient is specified for each stream- 
line in the subsonic region Equations (25) - (37) then yield flow proper- 
ties at all the grid points. This process is repeated for subsequent axial 
stations until Station II is reached. Station II is either a specified 
number of ax steps from Station I or is the station at which the Mach number 
along the upper bounding streamline falls below some prescribed value. 

At Station II, the mass flow contained between the upper and lower subsonic 
boundaries (GB) is evaluated as follows: 


I II 



FIGURE 13. 
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Referring to Figure (13), let RS denote the streamline below FG. Along RS 

2 


y Q = y R + (tane R ) (x $ -x R ) + (— |-) + f (x $ -x R ) 


cos"e 


X S -X R’ 


0 . 


tane s = tane R + ( — j-) (x $ -x R ) + j (x s -x R )‘ 

COS 6 


while between G and S 


1+J 1+J 

y _y 

+ \ ((pP c ose) s +(pqcose) G ) S 1 + j G — 


(39) 


(40) 


(41) 


where ^ ^ and = ip R 


Equations (39), (40) and (41) may be solved for y^, and e. This pro- 
cess is continued down to the wall (point B of Figure 12) yielding a 
* 

value of y B different than the value of y B obtained using Equation (18). 

A value d is obtained from Equation (18), which would adjust the lower 

★ 

wall shape so that y B =y B . The numerical procedure is restarted from Sta- 
tion I with a new value of d' chosen such that d‘ = . The iteration 

procedure for d appears to be self converging. Note that in repeating 
the calculation of the subsonic flow field ABFG with a new lower wall shape, 
the supersonic portion of the flow is not recalculated. Flow properties 
along FG are stored so that the derivatives can be made for evaluating the 
shear terms. It is assumed that the slight perturbations in the wall shape 
AB occurring in the iteration process do not influence the value of the 
terms stored along FG. Having converged on the lower wall shape, a new 
upper matching streamline. is determined .at Station II based on the Mach 
number profile. Then, the calculation proceeds from II to III as described 
above. This process is repeated until the sonic line closes and the entire 
flow field is again calculable by "viscous-characteristics". 
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VI. SHOCK PHENOMENA 


Numerical techniques have been developed for analyzing shock discontinuities 
occurring in a viscous combusting flow field. The analysis of higher order 
shock structure has been discussed extensively in References (16) and (17). 

For the problem under consideration, it appears that the modifications to 
the Rankine-Hugoniot relations to include the effects of heat conduction, 
viscosity and shock structure (which involve the local shock curvature) in- 
volve higher order approximations than the inclusion of these terms in the 
characteristic relations, hence the shock model employed uses the standard 
Rankine-Hugoniot relations, with the chemistry frozen across the shock. 

The flow may be nonuniform on both sides of the shock and "viscous- 
characteristics" are used in performing a shock point calculation, which 
takes the transport and chemistry terms into account. 

A. Shock Point Calculation - Assume a coordinate system oriented 
along (t direction) and normal to the shock surface (n direction) as shown 
in Figure (14). The angle beta ($) is the direction cosine of the shock 
with respect to the Cartesian direction x, and u and v are the velocity com- 

A A 

ponents in the n and t directions, 

(42) 

(43) 

(44) 

n = - sin fi f x + cos )3 { y 
t = cos /3 t x + sin /3 ty 
M oo = ' M n 4 + M t t 


FIGURE 14. 


n = - sin3i v + cosei 
x y 

t = cosei + sin&i 

X jr 

~ y * /x, « * 

V = - u n + v t t = ui x + viy 

Y 
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The Rankine Hugoniot relations take the form 

(45) 

(46) 

Tangential Momentum: v. = v. (47) 

z l z 2 

Energy: H = h + j = constant (48) 

where 

h = s cijh . (T) 
i=l 1 1 

State: p = p(P, T, a.) (49) 

The a. 1 s are the mass fractions (and remain constant across the shock) and 
h.j's the enthalpies of the i** 1 chemical species. Employing the jump rela- 
tions for a given shock angle and upstream conditions requires an iteration 
process since the mixture is calorically imperfect. 

Let 1 designate upstream conditions and 2 downstream conditions. To solve 
the jump relations knowing conditions at 1, a value for is assumed. The 
density using Equation (45) p 2 is computed; is computed using Equation 
(46) and Equation (48) yields a value for h 2 which can be inverted by a 
local iteration to find T 2 , since the composition is assumed frozen. The 
state Equation (49) then yields an alternate value for the density. If 
this value for density does not agree with that calculated from continuity 
to within a specified tolerance, a new value of is assumed and this pro- 
cess is repeated until convergence is achieved. 

Referring to Figure (15), a typical shock wave calculation is performed as 
follows. A value of the shock angle is assumed, which locates the point 
C. Since properties are nonuniform on the upstream side of the shock, a 


Continuity: 
Normal Momentum: 


0 . 


P 1 U 1 " P 2 U 2 


P 1 + p l u l P 2 + P 2 U 2 
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UPSTREAM OF 
SHOCK WAVE 



POINT 

POINT 

SHOCK 


FIGURE 15. 
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viscous characteristic calculation yields properties at Cl. The jump re-' 
lations (Equations 45 - 49) are solved using the determined upstream con- 
ditions based on the assumed angle This yields all . properties at C2. 
Using the deflection angle 8 q 2 calculated from the jump conditions, a 
viscous-characteristic calculation performed along (C2-A2) yields an al- 
ternate value of the pressure at P^ . The pressures are compared and if the 
difference exceeds a specified tolerance, a new value of is assumed and 
the process repeated until convergence is obtained. 


B. Under-Expansion Interaction - The program developed has the 
capability of computing the interaction produced by pressure mismatch between 
a jet and a surrounding airstream, for the case of an under-expanded jet. 
This situation is depicted in Figure (16). It is assumed that during the 
under-expansion interaction, the species remain chemically frozen. The ex- 
pansion is assumed to be isentropic and the local interaction is two dimen- 
sional and inviscid in the limit of vanishing radial distance with respect 
to the injector lip. 


The basic equations describing the Prandtl -Meyer expansion process are 


P/pT = constant 
1 2 

h + |- V = constant 
+ i d(V 2 ) = 0 

P <- 


- d 1n(P) ± d e = 0 

Y 

For a small incremental step aP, Equations (52, (53) and (50) can be 
written 


V o - V 2 = - ft y ) [ft ft ] 

L. 1 Y“i P 2 Pi 


(50) 

(51) 

(52) 

(53) 


( 54 ) 




FIGURE 


_ UNDEREXPANSION 
r SHOCK WAVE 

^-SLIP LINE P + =P~ 

e + =e~ 

PRANDTL - MEYER 
EXPANSION 
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1 in (Pg/pp ± (e 2 -e 1 ) = 0 

. (55) 

P 2 /p 2 = ’*Y P 1 

(56) 


where y is held constant in the integration step, yielding values for V 2 » 
p 2 » and e 2 . Then Equation (51) yields h 2 and an inversion yields T 2> In 
this manner, the Prandtl-Meyer equations may be integrated* 

Since the flow deflection and pressure downstream of the shock wave and 
Prandtl-Meyer are unknown an iteration process is required. A typical in- 
teraction calculation proceeds as follows. A shock wave angle is assumed 
for which flow properties (P,T,e etc.) are computed downstream of the shock 
wave. Equations (50) through (53) are solved using small increments of aP 
with the pressure behind the shock as the final pressure and the jet pres- 
sure as the initial pressure. If the turning angle for the expansion does 
not agree with the flow angle behind the shock to within a specified tol- 
erance, a new shock wave angle 3 is assumed and the process repeated until 
convergence is obtained. After this solution is obtained, the initial pro- 
file is adjusted by spreading the expansion over a small finite region and 
the program marches downstream carrying the shock wave as a discontinuity 
surface in the flow field. 

C. Shock Coalescence - The presence of significant compression 
waves generated either by combustion phenomena or physical boundaries as 
shown in Figure (17) requires that the entropy rise associated with the 
focusing of these v/aves be included in the calculation. Since the numerical 
scheme follows streamlines, not characteristics, a detection technique 
based on pressure profiles must be employed, rather than a crossing of waves 
of the same family. 

In strictly inviscid computations, detection of coalescence is based on 
determining the wave strength associated with the crossing of characteristics 
of the same family. However, the introduction of transport phenomena has 
the effect of dispersing the wave so as to weaken the usefulness of the above 
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RAPID COMBUSTION 
REGION 



SHOCK WAVE 


^-STREAMLINE 
DEFLECTION 
DUE TO COMBUSTION 



SCHEMATIC OF STRONG COMBUSTION WAVES 


FIGURE 17. 

detection criteria. A more practical approach is to locally determine the 
shape of the pressure curve once the above crossing has been detected. That 
is, a local polynomial of the form 

y = a + bp + cp 2 + dp^ (57) 

is fit using additional data from neighboring points in the region where cross- 

2 2 

ing is detected. Then, the simultaneous vanishing of dy/dp and d y/ dp indi- 
cates the presence of a coalesced shock wave. Typical pressure distributions 
in the region of strong waves are shown in Figure (18). 

The detection procedure is as follows. From data at the initial station the 
following equations for waves of the same family are solved for a crossing as 
shown in Figure (19). 


dx 


C, I 


= tan ( e ±y ) 


dx 


C, 1+1 


= tan (e±y) 


1+1 
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When a crossing is detected between points I and 1+1, a polynomial of the' 
form given in Equation (57) is fit using data at points 1-1, I, I+l, 1+2. 

If an inflection point exists for the above polynomial between I and I+l, 
it is determined if the magnitude of the slope y 1 at the inflection point is 
less than a specified tolerance. If so, an embedded shock wave is inserted 
between points I and I+l. The shock wave angle is assumed to be average of 
the characteristic angles, 

fi = ( (e±p) t + (e±p) I+1 )/2 (58) 

the + sign for an uprunning shock and the - sign for a downrunning shock. The 
inserted shock wave replaces the grid points I and I+l. The program allows 
for only one embedded shock of a given family, and cannot handle shock cross- 
ings or shock-boundary interactions. 



FIGURE 19. 
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VII. SAMPLE CALCULATIONS 

Case I : An axially symmetric hydrogen jet is injected into a Mach 2 air 

stream as depicted in Figure (20). 



FIGURE 20. 


Nominal conditions are 

Hydrogen - injection velocity 6850 ft/sec 
static temperature 190°K 
static pressure 2 atm 
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•Air - free stream velocity 4300 ft/sec 
static temperature 1255°K 
static pressure 1 atm 

The thickness of the injector lip and the effect of air side boundary layer 
leads to the formation of a base region downstream of the injector lip. Due 
to the large residence time in this region, combustion is assumed to initiate 
here. The pressure mismatch across this finite width region is accounted for 
by initiating the calculation a small distance downstream of the injector lip 
as depicted in Figure (21). The flow in the base region is assumed to be in 
chemical equilibrium and properties are obtained by assuming a smooth $ varia- 
tion between the hydrogen and air streams and obtaining the property variations 
by mass averaged formulas for the total enthalpy and velocity, with the pres- 
sure specified a priori. The program was run a total of 2000 axial stations 
(from x = 0 to x = 34). The overall flow field is depicted in Figure (22). 

The underexpansion shock is carried as a discrete discontinuity; its inter- 
actions with the upper constant pressure boundary (occurring at x = 10) was 
performed by a hand calculation. Representative flow streamlines, the iso- 
therms T = .5, T = 1.2 and T * 1.5, and the $ =1 line are also depicted in 
these figures. The wave field tends to significantly influence properties 
in the combustion zone as evidenced by the closing and reopening of the 
T = 1.5 isotherm. For example, this isotherm closes when the initial down- 
ward Prandtl -Meyer fan reflects off the axis and passes through the combus- 
tion zone (in the region 2 < x < 8) and again closes after the expansion from 
the upper surface crosses ('in the region 18 < x < 24). The variation of pres- 
sure, Mach number and mass fraction of hydrogen along the axis is depicted in 
Figure (23), A Ferri-Kleinstein eddy viscosity model was used for the flow 
in the potential core (defined by > .99) the viscosity varying from 
10" 4 lb-sec/ft^ at x = .5 to 2.14 * l!o~ 4 at x = 10. Downstream of this station, 
the value was maintained at 2.14 * 1G~ 4 . Flow properties are tabulated at 
the following stations: KOUNT =0; x = .5 

KOUNT = 100; x = 1.297 
KOUNT =500; x = 6 ,6346 
KOUNT = 1000; x = 14.932 
KOUNT = 1500; x = 24.97 
KOUNT = 1800, x = 30.24 
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CASE I, 

PROGRAM VIS - CHAR 
WITH 

EMBEDDED SUBSONIC FLOW 
SHOCK WAVES 
AND FINITE RATE H 2 - A I R CHEMISTRY 


TYPE OF FLOW IS “AX I SYMMETRIC 
CHEMISTRY IS FINITE RATE 
— JET"OR~NOZZLE" RADIUS' (RTH) - ; I3000E-01 FT • 


REFERENCE CONDITIONS 


MACH NO. .TEMINF) =” .19000E+01 

VELOCITY (UIN) = .43069E+04 FT/SEC 

TEMPERATURE" '(TIN)""- .12550E+04 DEGREES K 

PRESSURE (PRES) = .21160E*04 LB/FT**2 

DENSITY (RHOINF) = .54305E-03 SLUGS/FT**3 

FROZEN SPECIFIC HEAT RATIO (GAMINF) = .13199E+01 

MOLECULAR WEIGHT (WINF) ="" .28850E+02 ~ 

REYNOLDS NUMBER (RE) = .34450E+05 

PR'A HD T L”N U M B E R " " ( P R } “= * T 0 0 0 0 E + 0 1 *” 7 

LEWIS NUMBER (XLE ) = .10000E+01 


OUTPUT HEADINGS 


X - X/RTH 

y— — yy p T H 

0 - VELOCITY/UIN 

' T - TEMPERATURE/T I N — — 

P - PRESSURE/PRES 
TR-^te"OW" DEFLECT ION~< RADIANS) 

EM - MACH NUMBER 

GAM - SPECIFIC HEAT 

XMASS - NON-DIMENSIONAL MASS FLOW 

“ PHI -"EQUIVALENCE RATIO 

W - MOLECULAR WEI GHT 

MASS FRACTIONS 

~~ ALP ( V) - ' H “ “ 

ALP ( 2 ) - 0 

—ALP (3 ) - H20 — 

ALP (4 ) - H2 

~A'LPT5T~^~02~: 

ALP ( 6 ) - OH 
ALP (7) - N2 
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KOUNT = 

0 




X = 

* 50 QOOE^OO 





SHOCK TYPE 3 

BETA = • 662E^00 

VISCOSITY = .99380E 

-04 <LB«SEC/FT**2> 


PT. 

Y 

Q 

T 

P 

I 

• 6 0 0 0 0 E ♦ 0 0 

V15376EV01" 

* 1 51 50E + 0 0 

• 20 0 OOE ♦01" 

2 

• 694 0 OE ♦ 0 0 

• 15376E+01 

•15150E+00 

• 20000E+0 1 

3 

• 762 OOE + 0 0 

• 1 6489E ♦ 0 1~ 

~.13700E^0G 

• 14200E+01 ' 

4 

• 850 0 OE + 0 0 

• 1 2 1 24E ♦ 0 1 

•51400E+00 

.1^200E+01 

b 

#92000E+00 

•11 076E+Q1 

•95000E^OO 

• 14200E^0l 

6 

• 99 0 0 OE + 0 0 

• 1 0384E + 0 i 

• 15000E + Q 1 

•14200E^01 

— 7 

. 1 053 0E + OT 

'• 945 1 0 E ♦ 0 0 “ 

“ *21 OOOE + O 1 

•14200E^01 ' 

8 

• 111 00E+01 

. 92283E + 0 0 

. 2 1 600E* 0 1 

• 1420QE*01 

9 

* 1 1 8 0 0 E ♦ 0 1 

•94262E+00" 

• 1 97 0 OE ♦ 0 1 

• 1420 0E + 01 

10 

• 12500E+0 1 

.93400E+00 

•1460 OE +01 

• 14200E+01 

IT 

VI 3200E+0 1 

•91937E*00“ 

■ V 1 0850E ♦ 0 1 

• 1420 OE + 0 1 

12 

• 1 390 OE + 0 1 

•91937E^00 

•10850E+01 

« 14200E + 01 

T3 

“VI 3900E + 01 

• 10001EV01 

• 1 00 OOE + 01 

•10000E+01 

14 

•8250oe*oi 

•lQOQlE+Ol 

•lOQOOE+Ol 

• 1 OOOOE^Ol 


PT • 

ALP ( 1 ) 

ALP (2 ) 

ALP (3) 

ALP (4) 

r - 

Til OOOE-09 

V11000E-09 - 

V 11000 E- 09 - 

• 10000E + 01 

2 

•11000E-09 

• 1 1 000E-09 

•11000E-09 

« lOOOOE^Ol 

V 

• 1 1 00 0E-O9 

• 1 1 OOOE-09 

•11000E-09 

• 1 0000E + 01 

4 

• 1 1000E-09 

•11000E-09 

•16300E*00 

•34000E^00 

3 

•T1000E-09 

•I 1 OOOE-09 

•POBOOE^OO 

• 18400E*00“ 

6 

• 1 1 OOOE-09 

• 1 1 OOOE-09 

•23300E^OO 

•7800QE-01 

7 

el 1 00 0E-09 

• 11000 E- 09 ~ " 

. 23700E* 0 0 

V13000E-01 

8 - 

• 1 1 000E-09 

•24000E-02 

•21200E+00 

•11000E-09 

9 — 

•11000E-09 

•20000E-02 - 

* 1 450 0 E + 0 0 

• 1 1 000E-09 

10 

• 1 1 OOOE-09 

.1 1000E-09 

•52000E-0 1 

•11000E-09 

n — 

TITO 0 0 E- 0 9 — 

11000E-09 

.1 1 0 0 0 E - 0 9 

•110G0E-09 

12 

• 1 1 00 0E-09 

.1 1000E-09 

•11000E-09 

• 1 1 OOOE-09 

IT 

a 1 000E-09 

V1 1000E-09" - 

•11000 E- 09 

• 1 1 0 00E-09 

14 

• 1 1 OOOE-09 

• 1 1 OOOE-09 

• UOOOE-09 

• 1 1 OOOE-09 
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TH 

EM 

RHO 

•92249E+00 “ 

•92249E+00 

•72429E+00 

• 48997E + 0 0 
•41609E+00 
•43Q52E ♦ 00 
•5Q4Q2E+00 
•56660E+00 
. 64958E+00 

• 93850E + 00 
• 13086E*01 
• 13086E+Q1 
•99992E+00 
.99992E+00 

GAM 

XMASS 

5 cci i c ♦ n ft 

0 . 

0. 

•10450E+00 ‘ 

• 1 0450E+00 
.10450E+00 
.10450E+00 

• 1 0450E + 0 0 

• 10450E+00 

• 1045GE+00 

• 1 0450E + 0 0 
. 1Q450E + 00 

• 10450E* 00 

o, - - 

0. 

• 19100E+Q1 ' 
•19100E+G1 

* *21500E + 0l 

• 13200E+01 

• 1 1300E+Q1 
•11000E+01 

•iioooe*oi 

• 1 1 40 QE * 0 1 
•12400E+01 
•14600E+01 

• 1 680 OE + Q 1 
. 16800E + 01 

• 19000E + 01 

•19000E*01 

• I 424 l 
• 14243E*01 
• 14295E+0 1 

• 1 3869E+0 l 
•134J4E+Q1 

• 1 2873E + 0 1 
• 1 2485E+0 1 

• 1 2459E* 0 1 

• 1 2596E *01 
• 1 2688E*0 1 

• 13151E+01 
• 13151E+01 
• 1 3200E+0 1 
• 1 3200E+01 

•34158E+00 
-.-4O608E^OO— 
•469 16E+00 

• 50 165E*00 ■- 

•53183E+00 

•56138E+00 — 

•592Q1E+00 

•63726E+00 - 

•70023E+00 
•79325E+00 - 

•90674E*00 

- • 90674E + 0 0 

• 3397 1 E*02 



ALP < 5 > 

ALP (6 ) 

ALP (7) 

PHI 

W 

• 1 1 000E-09 

•11000E-09 

-• 55000E-09 

-.26620E+12 

• 20 160E + 0 1 ~ : 

• 1 1 OOOE-09 

• 1 1 00 GE“09 

-•55000E-09 

-.26620E*12 

•20160E+01 

•11000E-09 

•11000E-09 

-.55000E-09 

-.26620E+12 

• 20 1 60E + Q 1 

• 1 1 000E-09 

.11000E-09 

• 4970 OE ♦ 0 0 

• 19142E+02 

.511 67E*0 1 

• 1 1 000E-09 ' 

.11000E-09 " 

* 60800E + 00 - 

— *89665E*0 1 

- • 80309E* 0 1 

• 1 1 00 0E-09 

•11000E-09 

•68900E*00 

• 39 8 3 5E* 0.1 

• 1 3 1 20E + G2 

• 1 1 00 OE-09 

•55000E-02 

• 74450E* 0 0 

^14231E*01)" 

•21504E+02 - - 

.15000E-01 

• 14000E-0 1 

• 75660E + 00 

a 863 1 6E + 0 0 

•24865E+02 

•aooooE-oi 

•85000E-02 

. 7645QE^00 

•58345E+00 • 

•25999E+02 

•181ooe*oo 

•11000E-09 

•76700E+00 

•2007 1E + 0Q 

* 27838E*02 

. •232G0E+00 

•11000E-09 

• 76800E«-00 

•81900E-08 

• 28848E + 02 

.23200E+00 

•11000E-09 

• 76800E^0 0 

• 61 90 QE-08 

•28848E+02 

.23200E+00 

•11000E-09 

• 76800E + 00 

• 8190 0E-08 

- 288486*02 

• 232G0E + 0 0 

• U0OOE-O9 

•76800E+00 

•81900E-08 

•28848E+02 


TRT69 

KOUNT = 100 


Page bO 


X = 


.12970E*01 


SHOCK TYPE 3 BETA a 


• 645E+00 


VISCOSITY = .11 356E-0 3 (LB»SE C/FT »»2) 


PT. 

j- 

2 

3” 

4 

5 _ 

6 

7~ 

8_ 

9 

10 

IT" 

12 

rr 

14 

— 15 - 

16 

IT" 

18 

rr 

20 

2T 

22 

23" 

24 


Y Q T 

— 0“ . 1 5’8 0 2 E + OX ' V 1 4 4 8 9 E ♦ 0 0 

• 1 0326E + 0 0 * 1 5839E + 0 1 .14432E+00 

. 20753E + 00" • 1 59<+3E ♦ 0 1 .14267E+00 

• 31 376E + 00 . 16099E + 01 .14017E+00 

.42282E+00 .16284E+01 .13718E+00 

.53535E+00 • 1 6463E ♦ 0 1 .13428E+00 

VSS 15 OE VDTT . 16587 E ♦ 01 . 1 3 2 7 4E ♦ 0 0 

• 76496E + 00 .16526E+01 .14005E+00 

. 85728E + 0 0 .15847E + 01 il8896E+00 

. 95107E+00 • 1 347 7E * 0 1 .37483E+00 

- . 1 0 1 02E + 0 1 .12266E + 01 .62508E+00 

. 1 0699E + 0 1 .1 1 345E ♦ 0 1 .98895E+00 

V 1T273E + 01 VI 0692E + 0T~~~~* 1 3772E+0 1 

.11 820E + 0 1 .10257E+01 .16631E+01 

.12490E+01 .99254E+00” .17183E+01 

• 131 94E + 0 1 .97145E + 00 .14802E+01 

# X3956E^01 . 95295E + 0 0 ' .11800E+01 

.15023E+01 • 94849E + 00 .10818E + Q1 

7160 3 8 E + 01 79 4709E+0 0 VI 0 6 8 2 E ♦ 0 1 

• 1 7052E + 0 1 .94499E+00 .10658E+01 

. 18Q79E + 01 .94176E + 00 . 1 0 671E + 0 1 

.19941E+01 .93222E+00 .10753E + 01 

71 9941 E ♦O 1 • 10O01E*0l ‘ .10000E+01 

.82500E+01 .10001E+01 .10000E+01 


P 

•17290E+01 
. 1 7069E ♦ 0 1 

• 1 6434E + 01 
. 1 5499E* 0 1 
. 14428E ♦ 0 1 
.13421E*01 
.12689E+01 

• 1 2477E *0 1 
.12600 E* 01 
.12359E+01 
. 1 2348E+0 1 
. 12403E ♦ 0 1 
.12483E+01 
.12554E+01 

• 1 2622E+0 1 
. 12695E + 01 

• 1 2785E+0 1 
. 1 2643E+0 1 
.12679E+01 “ 
o 12779E+01 

• 1 295 1E + 0 1 
.13473E*01 
.lOQOOE+Ol 
.10000E+Q1 


PT. 

— 1— 

2 

— 3 _ 

4 

"5“ 

6 

r _ 

8 

— 9 
10 
IT" 
12 

14 

“15“ 

16 

”17- 

18 

20 

“ 21 " 

22 

23 


ALPU) 

TTTOO0E= r O9 

*110 0 OE-09 
V1 1000E-09- “ 
.U016E-09 
". 1 1671E-0 9 ' 

• 36806E-0 9 
T90935E-08 

.25565E-06 
“.43251E-05 
. 480 09E-04 
”• 1 9326E-03 
.5591 IE-03 

-T12743E--02 

. 20499E-02 
10425E-02” 
.15203E-03 * 
—.99 1 97E-05 
.31768E-06 
“V28359E-07 
• 94356E-08 
“.46227E-08 “ 
.11000E-09 

• 1 1 0 00E-09 


ALP ( 2 ) 

-VlTOOOE-09" 

• 1 1 OOOE-09 
"V11O00E-09- 

.11000E-09 
“VI 1000E-09 
.11011E-09 
“.IT 327E-0 9 ” 
.18863E-09 
” ;T2752E-08- 
.1 17 78E-07 
“.-725S6E-07- 
.18134E-05 
“.49932E-04 
.87296E-03 
-.61256E-02 
.31695E-02 
”.42284E-03 
.25394E-04 
V32632E-05 
. 12869E-05 

• 62972E-06' 
. 1 1 QQQE-09 
"• 1 1000E-09 


ALP t 3 ) 

• 1 1353E-09 
. 1 534 3E-0 9 
© 1 31 75E-08 
.35177E-07 
.93569E-06 
.21702E-04 

.41 574E-03 
. 56903E-02 

• 3849 1 E-Q 1 

• 1 3375E+00 
. 17892E + 00 
. 206148*00 

.21 689E ♦ 0 0 
.21 146E+00 
» 15255E+0G 

• 72325E-01 
.1.7624E-01 
.326G2E-02 

' .12109E-02 
.57934E-03 
- .28608E-03 

• 1 1 000E-09 

.11000E-09 
i i n n n r * ri r» 


ALP ( 4 ) 
•10000E+01 ■ 
• 10000E+01 
. lOOOOE^Ol 

• 1 OOOOE^O 1 
.10000E+01 
.99991 E ♦ 0 0 
.99833E+00 

• 977 1 2E^0 0 

• 84542E + 00 

• 46248E ->, 00 
. 27374E ♦ 00 
. 14340E^OO 
.65675E-01 
. 22949 E -01 
. 35832E-02 
•44560E-03 
.45960E-04 
.26649E-05 
.25785E-06 

• 8581 3E-07 
.41584E-07 
. 1 1000E-09 
.11 OOOE-09 

\ innn r o n 


a r 
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\ 


TH 

EH 

RHO 

GAM 

XHASS 

0. 

•20056E+01 

".83388E + 00 

• 14266E ♦ 0 1 

0. 

•16342E-01 

•20141E^Q1 

* 82644E^00 

ol4268E*01 

•70014E-02 

"V33855E-01 

“.20385E*01 

. 80488E*00 

• 14274E^0 1 

- v27995E-01 — 

•52889E-01 

♦ 20762E + 0 1 

•77267E*00 

• 14283E+0 1 

• 62945E-0 1 

.72650E-01 

•2122QE+01 

•73499E^00 

* 14294E^01 

• 11 187E*Q0 

. 90830E-0 1 

•21676E+01 

•69849E+00 

« 14306E^01 

•17491E+00 

" .10434E + 00 

•21977E*01 

*669QOE^OO 

•14312E*01 

*25241E*00 — 

•11324E>00 

•21550E*01 

•63590E+00 

• 14282E>0 1 

•3387 lE^OO 

• 1 2140E^O0 

• 191Q2E+01 

.5430 1 E*QO 

* 14 1 32E^0 1 

• 40979E + 00 

• 1 1988E + G0 

• 151 1 8E + 0 1 

*45483E^0Q 

• 13935E+01 

•47180E*00 

•10812E*00 

• 13206E*0 1 

e 41393E^OO 

• 13778E^0 1 

•5041 OE+OO - 

.94879E-01 

• 1 231 OE^O 1 

.40895E^00 

•13345E^01 

•53412E>00 

' .85926E-01 

- .12176E^01 

• 44002E* 00 

• 12950E^01 

.56346E + 00- 

*8 1656E-0 1 

• 1 2509E+0 1 

•49906E^00 

• 12736E^0 1 

•59433E+00 

— »80118E“01 

— *1 3404E+0 i 

- .61396E+00 

- • 12709E + 0 1 

- •63984E+00 - 

•81569E-01 

• 14883E^0 1 

• 80428E^Q0 

• 12862E + 01 

•70253E^00 

— .87625E-01 

• 1 6634E* 0 1 

•'10685E + 01 - 

— • 13070E*0i 

- •79521E*00 — 

•93582E-01 

• 1 7340E^0 1 

• 1 1658E^0 1 

• 13146E*0 1 

•95877E+00 

•93206E-0 1 

• 1 7430E+ 0 1 

• 1 1 858E + 0 1 

;13158E+0r 

• M337E + 01 

.93062E-01 

• 1741 3E+01 

• U984E*01 

• 13160E+01 

• 1 3220E*0 1 

“ .94949E“01 

• 17345E+01 

• 1 21 33E+0 l 

• 1 3160E + 0 1 

•15263E+01 

•10422E^00 

•17108E+Q1 

• 12529E+Q 1 

♦13156E+01 

• 1 9332E + 0 1 

- .56521E-12 

•19000E^01 

• 99992E^00 

• 13200E^0 1 

• 19332E^01 

0. 

•19QQ0E*01 

•99992E>Q0 

• 13200E+01 

•33975E*02 


/ 







ALP (5 ) 

ALP (6) 

ALP ( 7 ) 

PHI 

W 

• 1 1 000E”09 

.1 1000E-09 

-.53926E-09 

-♦29833E* 1 2 

• 20160E^0 1 - 

•11000E-09 

•11000E-09 

-•41810E-09 

.82356E+12 

•20 1 60E^0 1 

— .11002E-09 

- .1 1000E“09 

• • 31 156E-08 

• 74400 E> 10 

.201 60E*0 1 - 

• U110E-09 

•11000E-09 

•10585E-06 

• 24956E+09 

• 20 160E*0 1 

•16051E“09 

•11005E-09 

- • 2836 IE-05 

• 93505E*07 

•20 160E*01 — 

•22714E-08 

•11128E-09 

•65739E-04 

• 40333E + 06 

•20162E+Q1 

— .84238E-07 

— •14003E-09 

• 12578E-02 

.2104 1 E*05 

*20191E^01 — 

•26809E-05 

•64085E-09 

•17187E-01 

. 15074E + 04 

•20593E*01 

- .50438E-04 

•50984E-08 

• 1 1603E+0Q 

• 1 939 1E + 03 

•23494E^01 

•60521E-03 

•25731E-07 

•40312E+00 

• 31 338E + 02 

•39796E^01 

•24574E-02 

•51979E-06 

*54470E*00 

• 14277E + 02 

•60450E*01 - 

• 65669E-02 

•22442E-04 

• 64332E+00 

• 68762E + 0 1 

• 94076E + 0 1 

-.12758E-01 

.47011E-03 

- • 70288E+00 

.34433E+01 

« 14006E*02 - 

•22337E-01 

•40222E-02 

•73631E*00 

•17631E+01 

.19074E*02 

•69350E-01 

•10785E-01 

. 75656E+QQ 

•78345E*00 

•24114E^02 

• 15522E+00 

•39004E-02 

• 76479E'*- 00 

.3Q871E^00 

•27053E*02 

•21 394E*00 

•50789E-03 

•76745E^00 

• 70723E“01 

• 28451E «*02 

•22874E+00 

•50370E-04 

• 76792E^00 

• 1 2720E-01 

*28779E«-02 

'"•23080E + 00 

•12445E-04 

- •76797E^00 

•46822E“02 

” •28823E^02- 

• 23143E+00 

• 56240E-05 

•76799E^00 

•22380E-02 

• 28836E^02 

•23172E^00 

• 2 7697 E “05 

• 76799E^00 

• 1 1050E-02 

• 28842E^ 02 

•23200E+00 

• 1 1 000E“09 

. 768GOE+QO 

•81900E-08 

•28848E^02 

.23200E+00 

•11000E-09 

•76B00E+00 

•81900E“08 

•28848E ♦02 

• 23200E+00 

•11000E-09 

.76800E*00 

•81900E“08 

•28848E+02 
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KOUNT = 500 


20 

-21- 

22 


X = 

• 66346E + 01 


i 



SHOCK TYPE 3 

BETA = .630E+00 

‘ VISCOSITY = , 1 7496E 

-03 (LB*SEC/F T **2 ) 


PT. 

Y 

0 

T 

P 

1 

2 

o. 

.23219E+00 

• lb/cHE + UT 
. 16777E+01 

• 1 2 94 0 1 + oO 
. 12999E+00 

• 1 1 64 4 1 +01 
. ) 1 60 3E +■ 0 1 

3 — 

4 

.469 I 9E + OD 

. 72282E + 00 

71^T697EWr 
. 16160E+01 

7TT4S6>_ +00 
. 16209E+00 

. I T4 7 6E+'Ui 
. 1 1269E + 01 

“ 5 

6 

796 1 8 7 L + 00 " ' 

. 11063 E +01 

•13139E+01 

.2591 or +00 
.39317E+00 

.10 97 otT U"1 
. 1 0850E+01 

r~ 

8 


• 1 23636+01 
. 1 1466E+01 

.52^+0 6T5 +"crcr 
.77879E+00 

• i U 7 8ttt + U i 
, 1 0729E + 01 

9 

10 

. 140 /6E + 01 
• 1 5353E + 0 1 

VTO 9S3ET+""0T' 
. 1 04S4E + 0 1 

.10 055E + 01 
. 13483E + 01 

•10o9bt+01 
.1067 6E +01 

12 

7TTTT5 1 E + 0 1 
. 18287E+01 

.TO 0"82E7 r 0T 
. 99714E+00 

• ljth jt + ij i 
. 11295E + 01 

• 10 fUlL + Ul 
. 10521E + 01 

T3 — 

14 

.215876 + 01 

• 2584 RF + 0 1 

.9991 1E+0U 
. 10002E+01 

1 r ) 1 0 5 E + 0 1 

. 10003E + 01 

• 9930 2E + 0 0 

1 5 
16 

. 30 \35T^~0-]- 

• 3 4 4 5 4 E + 01 

v 998‘6'4'E _ +~fr(r 
• 9945 1 E + 00 

.1 01) ID L + D 
.10064E+01 

• 1 Utl c CZ. + O' Ji 

•10225 E *01 

T7 

ie 

.39 723E+01 
• 44955E + 0 1 

. 98 8 04 1 + 00 
.97997E+00 

. 1 0 l J7E + 01 
. 10227 E + 01 

• 10 55 7E + 01 
. 1 096 /E + 01 

T9 — 

•5nuV5t+ol 

.96449 ET+TTO" 

• 10 J44E + 01 

. 1 I 57T6 E + tn 


•59174&+01 

T^rtTtr^rfjrt 
• 8250 OE + 0 1 


.94327E+G0 .10634E+01 .12660E+01 

TiDtre-i-E ^ 0 1 HOOTO E-+1H HOD O 0E+ OH 

• lOOOlE + Ol HOOOOE + Ol HOOOOE + Ol 


PT. 

— 1- 
2 

3“ 

4 

— 5“ 
6 

8 

10 

Tf 

12 

T3" 

14 

T5~ 

16 

TT~ 

18 

-TD- 

20 

- 21 “ 

22 


ALP(l) 


• 30 935E-05 
24^7 6f~- t)^r 
. 19358E-03 
T9T0T9E-1 


ALP (2) 


~=DT 
•45433E-09 
tTT 
•23481E-07 


ALP (3) 
t22T1TE'-D3~ 


• 8 1 657 E - 03 

•29465E-01 
^DT 


ALP (4) 

• 999076+" 00 — 

• 9966^+E + 00 . 
r97-tf filT^D — 

• 87648E+00 . 
v58-&84fr*DD — 


•18505E-0? , 78039E-07 H^2B9E+00 •38524E+00 

. 26118 6^02 r3 4~76D E ^06 H68T2 fr+TD •273888+00 — 

•38366E-02 .61792E-G5 .19656E+00 .15182E+00 

.4B409E-02 r * 83 2 ? E * Q 4 H H 3 38tTD 0 ra5r(r 7 76 - 01 

•56570E-02 .23471E-02 .21968E+00 .242808-01 • 

H 5 n 3 66^iH HD7 88E"+1r0 r?T6-59 E— -02 — 

•33364E-03 .81S82E-02 .32929E-01 .21706E-03 

.4226 0 E - 03 H-9 -H t 8 6- 92 r9 43516-0-5-— 

•30578E-06 .110556*04 .603636-04 .22668E-06 

,551905*06 *220356-06 — : — H95j42fr-05 rDB-frD I L " - 06 — 

• 18005E-0 9 *3214 7E-08 .244795*07 . 17070E-09 

Vl 1 T T2E -^ D D r42V59e-09 HD 0606-09 — 

•11004E*09 H 1204E-09 H5277E-09 .11005E-09 

“Til ODTfTtrO . 1 1 O94E-09 • 1 J<>6i 6-09 r HO 02 £— 09 — 

•11000 E- 09 H 1000E-09 .lloOOE-09 . 1 1000E-09 

T l ' l O ' CO E ~ 0 1 • 1 I 0 U 0 E - 0 9 * 1 i OOOE-09 H-tDD-fr E - 0 9 — 

•1100 OE -09 .11000E-09 H 1000E-09 H 1000E-09 
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.20237E-02 
-73861536^02- 
.661 00E-02 
78224 06-02“ 
.682926-02 
T5-2-9(FfE>02- 
.34106E-02 
7338 9 OE^O 2” 
.581376-02 
Tt5&03Mt- 
.21484E-01 
"Tt^5i>6E"0 1 
.138466-01 
-r34434E~03- 
. 1 9204E-0 1 

72 r& 826- or 

. 3837SE -0 1 
-rF26~3FE- 0F 
.88073E-01 
-*-43-743 6—32 


EM 

• 224T6T4'0T 

• 2247 16 + 01 
-.F2183E + 01 
. 2 0636E + 0 1 
-rr7547E + at- 
. 1 5465E + 0 1 
73^3376*33 
*135116+01 
7333026+03 
. 13489E+01 
i“l'5 537 E'+ 0 i - 

• 1 751 OE + 0 1 

• 1 89986 + 0 1 
-*-l-8-964e+01- 

.188386+01 
-rr865iE>-ar 
. 18420E+01 
Tt 8324£ _ +‘03 
.174036+01 
-.-390006+01- 
.190006+01 


, 625676+00 
T6O8O9E+00- 
,547766+00 
7 4 74 3 8 6+00““ 
.43958E+0Q 
742 51376+03” 
.425556+00 
73407 9 E + 0 0 
.483846+00 
-.6902 1 E-+ 00— 
. 892296+00 
-*^93536+00- 
. 99253E + 00 
ilOOHE+Ol- 
. 1 0 1 60E + 0 1 
73 04l“4E _ +0i _ 
.107236+01 
“1332 2 E + 0 1 - 
.12110E+01 
i-999926+OO- 
.999926+00 


GAM 
. 14325 

• 14323E+0 1 
V14-303E + 01' 

• 142046+0 1 
-rl-40 1 4E ♦ 0 1- 

• 1 39356+01 
73-38706+03 
.136126+01 
-733316E+03 

• 12970E+01 
-.-1-297 5E+ 01- 

. 13106E+Q1 
-*333926+03 

• 131996+01 
1-1--3199E + 01 

• 1 31966+01 

-.13T91E + 01 
.13I86E+01 
-.-1-31796+06 
. 1 3 163E + 0 1 
-.132006+01- 

• 13200E + 0 1 


•283736-01 
733 4396 + 00— 
• 25783E ♦ 00 
7-436086+00 — 

• 5 1 047E + 00 
7567436+00 — 

• 6448 1 E + 00 

7 704516+00 — 
.79738E+00 
-.962716+00 — 
♦11357E+01 
-*335426+03 — 
.27562E+01 
-*-395176 + 01 — 
.535346+01 

.7345 66+0 1 

. 96477E+0 1 
1-122356+02 

• 1 7738E + 02 
-.-677386+02 — - 
• 34262E+02 



ALP ( 5 ) 
-*-45928 E~ 05- 
.192036-04 
-rt624-06-03- 
.113516-02 
v5 2 133 E— "02- 
. 94 74 1 E-02 
-.319396-0-1- 
. 14013E-01 
-*3-3-7996—03- 
• 1 6338E-0 1 
-*3 0 943 E + 0 0- 
. 19125E+00 
-*-229796+00- 
.231946+00 


ALP ( 6 ) 
i339846— 09- 
.23304E-09 
“1-0-6546—0-8- 
.536736-08 
^921336-08- 
.292346-07 
-i-5 12476—06- 
.20472E-04 

.425656-02 
-s-7-4650 6—02- 
• 20289E-02 
*3-04296-03- 
.289966-05 


ALP (7) 
-696946-03- 
251936-02 
6627 6 E — 03 — 
92705E-0 1 
-30 8366+00- 
46054E + 00 
r5 43446+0 0 — 
, 63372E + 00 


PHI 

-.-379296+05- 
• 1 0473E + 05 
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Case II : An axially symmetric Mach 1.2 hydrogen jet is injected into a . 

Mach 3 airstream as depicted in Figure (24). The initial conditions for 
this case are tabulated under KOUNT = 0. The pressures of the jet and 
airstream are balanced and the effect of the initial airside boundary 
layer and base region are analyzed by assuming chemical equilibrium to 
prevail in this region as described for Case I. The combustion generates 
a mild compression field, the pressure increasing to about 1.25 of the 
initial value in a distance of about 2 slot heights. The axial pressure 
distribution for several streamlines is depicted in Figure (25). Note that 
initially the pressure is highest in the combustion zone, but when the com- 
pression waves generated by combustion reflect off the centerbody, the 
pressure becomes highest at the lower boundary. The overall flow field 
properties are tabulated. . The flow becomes subsonic at the axial sta- 
tion X = 2. Details of the sonic line and streamlines in this region are 
depicted in Figure ,(26).. Note that the combustion zone is in the vicinity 
y = 6, hence the streamlines are closer together in this region so that 
the mixing can be accurately calculated. 

In the region 2 < X < 2 * 2, the upper subsonic streamline (<j> = .591) is 
specified by the coefficients A thru D specified by the values of y, e, e s 
and e s$ at the station KOUNT = 365, while downstream of X = 2.2 a higher 
order term is added to turn this streamline down and hence re-accelerate the 
flow. The Ferri-Kleinstein viscosity model is also employed for this case, 
the viscosity varying from .6643 * 10 ^ lb.sec/ft^ at X - 0 to .8252 * 10 ^ 
at X = 2.4. The flow field properties are tabulated at the following stations 
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The flow field is calculated by “viscous-characteristics" from station 
KOUNT a 0 to KOUNT = 365, while the subsonic routine is used to calculate 
the flow between KOUNT = 365 and 405, below streamline 11. The flow down- 
stream of KOUNT = 405, is calculable by "viscous-characteristics" provided 
that the lower wall continues slope downward, hence, accelerating the flow . 
Note that the sonic line does not reach the wall since at station K0UNT=365, 
the wall is turned down hence accelerating the streamtubes in the vicinity 
of the wall. The wall pressure (as depicted in Figure 25) in the subsonic 
region continuously decreases as the wall turns down., consistent with the 
flow being supersonic adjacent to the wall. Hence, the neccesity of pro- 
viding a variable pressure in the subsonic region, employing the normal 
momentum equation is quite evident. A constant pressure scheme certainly 
would not provide the details of the flow field as obtained in this analysis 
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VIII. CONCLUSIONS 

The ’Viscous-characteristics" method for the analysis of supersonic viscous 
combusting flow fields has been extended to analyze local embedded subsonic 
regions as well as to analyze shocks produced by combustion and/or pressure 
mismatch between the injected gas and air stream. The numerical technique 
developed for analyzing subsonic zones requires that the subsonic region be 
bounded by a lower wall and an upper boundary that is slightly supersonic. 
The shape of the upper boundary is fixed by the flow conditions by requir- 
ing that y^e,. and e be continuous, but higher order terms may be sped- 
fied to shape the streamline arbitrarily. The shape of the lower wall is 
a function of the upper wall shape specified. The program developed can 
be used in its current form to design centerbody shapes for a specified 
streamline shape in the combustion zone. Note that the program can be 
modified to make the upper boundary a specified pressure boundary and hence 
a lower wall shape can be obtained fora specified pressure distribution in 
the combustion zone. The program can also be extended to analyze embedded 
subsonic zones surrounded by supersonic flow on both sides provided that the 
flow can be assumed inviscid on one side. Referring to Figure (7), where 
the flow beneath the embedded region is inviscid, a marching scheme can be 
developed wherein the flow in the subsonic zone and above is analyzed by 
a mixing type grid, while C + characteristics are followed in the inviscid 
region. It is felt that these program modifications would be realistic 
extensions of the current effort and should be considered in plans for 
future research in this area. The extension of this program to embedded 
zones with viscous supersonic regions on both sides appears to require a 
more significant effort. 

While the "viscous-characteristic" program developed for supersonic flow 
fields can be readily run by one with only a limited background in gas- 
dynamics, it is felt that the subsonic version be run by an analyst who is 
thoroughly familiar with the physics of the problem as well as the details 
of the analysis included in this report. 
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APPENDIX I 


MARCHING SCHEME - FINITE DIFFERENCE METHOD 



FIGURE 1-1. 

Referring to Figure (1-1), all properties are known at the initial station I 
and hence all derivatives 3/ay may be calculated at the mesh points 1, 3, 5 
etc. The derivatives p and e are related to the streamwise derivatives p s 

y y 

and 0 by the relations 

2 2 2 

cos 6 [A] yPM cos 00“ sine p 

p = ? o ^ 

s frcos^e - 1 


(AI-1) 
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e s = 


2 * 2 
sinecose [A] - (M -1) cose P - sineyPM ey 

JL — — 

yPM 2 (M 2 COS 2 0 - 1) 


(AI-2) 


where 



Jsine + 
y yPM 2 


(y-i)s 2 

ypq 




(y-1) h i * i 
(y -l)M 2 yPq 

V CO ' CO 9 * 


w 

pq 


£ 



(AI- 3 ) 


Assuming a value of e (or prescribing this value as a function of x in 

So 

the inverse problem), all properties at 2 may be obtained by the boundary 
characteristic relation along 2-B and the s-momentum, energy and species 
diffusion equations applied along the streamline 1-2, 


Knowing 
and (2). 2 


and 0 , 

s 2 


we may obtain P 

y 2 


and 0 


by inverting Equations (1) 


Combining the s-momentum, energy and species continuity equation, the rela- 
tion 


(H 2 -l) P s + yPM 2 0 n = [A] 

is obtained, while the normal momentum equation is written 

P + yPM 2 6 = 0 
n s 

Between points (2) and(4) 


P 


4 



+ [ \ +p y ] 


Ay 24 

2 


0 4 = 0 2 + 


[0 +0 ] 
y 4 y 2 


Ay 


24 


(AI- 4 ) 


(AI-5) 


(AI-6) 


2 


(A I -7 ) 
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Between points (3) and (4) 


AS 


P 4 = P 3 + C p s 3 +P s 4 3 -Z 


34 


(AI-8) 


AS 


Q = 9 + [e +e ] 
43 s 3 S 4 c 


34 


(AI-9) 


The system of Equations (4) to (9) (with Equations 4 and 5 applied at point 
(4)) along with the transformations 


and 


Is = cose + sine ty 


|_ = cose |y - sine fj- 


can be combined to yield the relations 


2 ( 64 - 63 ) , z 2 


4 * \ - y 4 M 2 S3 (M 4 2 cos 2 e 4 -X) 


2 Y 4 M 4 (© 4 “ 02 ^ s ' 1 * n6 4 


Ay 


24 


Y 4 M 4 0 y 2 s1n0 4 + 


2(MJ-l)cose 4 


\ 


Ay 


24 


(AI-10) 


2(M4“l)P 2 cose 4 


a ^24 " ^ 4_1 ^ P y 2 cos9 4 ” Sine 4 cos6 4 [A] = 0 


and 


? 2 , ^/VV- 0 ?) 0050 /! ? 

P 4 * j 2(H 2 cos 2 e 4 -l) ♦ i - T4 M 4 COS0 4 y 2 


2sine>, \ 2P~ 9 ? 2P„sine 

+ 2- I - (P + -rM (MfcosS,-!) - 


Ay 24 


s 3 &s 23 ' "'4^ 4 


Ay 24 


(AI-11) 


- P 


' sin8 4 - [A] cos 2 e 4 


= 0 
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Equations (AI-10) and (A I - 1 1 ) both are of the form 


P*f(e,M) + g(e,M, [A] ) = 0 


(AI-12) 


An iterative procedure for the solution of the above system proceeds as 
follows : 


(a) 

(b) 

(c) 


A value for P 4 is assumed (i.e., P 4 = p 3 + p as,.). 

3 

Application of the s-momentum, energy and species equa- 
tion along 3-4 -+■ q^, T^, a. M^. 

Equation (AI-10) -> 9^ by an iterative procedure. 


(d) Equation (AI-11) -> P^ which is compared to the value of 
P^ assumed. 

•k 

If I P 4 ” P 4 I > where z is some specified tolerance, a new value of P^ is 
assumed and steps (b), (c) and (d) are repeated until convergence is obtained. 


The calculation y points 6 , 8 and 10 are performed in an analogous fashion, until 
the lower boundary 9-10 is reached. In the direct problem 9-10 may be a wall, 
axis or lower matching streamline (in which case a P-e relation exists along 
the C + characteristic 10-A), while in the inverse problem 9-10 must be a wall. 

The direct problem requires that the deflection angle e 1Q obtained by the sub- 
sonic solution match the wall angle 8 W , = 0 if 9-10 is an axis, or satisfy 
the compatibility relation along A-10. If it does not, the value of e s must 
be iterated upon and the entire system solved repeatedly until the value of e s 
that satisfies the lower boundary conditions is obtained. In the inverse 
problem, the value of e 1Q determined by the subsonic solution yields the shape 
of the lower boundary and no iteration is required. 
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APPENDIX II 

MARCHING SCHEME - POLYNOMIAL METHOD 

In the subsonic region, the pressure distribution and flow inclination are 
expressed by the following power series in y: 

p(x,y) = Pq(x) + Pj(x) y + P 2 (x) y 2 + P 3 (x) y 3 * 

e(x,y) = e 0 (x) + e x (x) y + e 2 (x) y 2 * 

Referring to Figure (II-l), the subsonic region is bounded by the two stream- 
lines DjC^and D^C^ along which the Mach number has low supersonic values. 

The subsonic region is matched to the supersonic regions employing the 
viscous characteristic compatibility relations along the characteristics 
BC 2 and ACj. The analysis developed also considers the possibility that the 
subsonic region extends to an axis or lower wall. 


(AII-1) 

(AII-2) 


The derivatives P s and e $ at points and are calculated employing Equa- 
tions (AI-1) and (AI-2). The following procedure is used to calculate the 
subsonic region. 

Step 1. A value for is assumed (or prescribed for inverse 

Sr 

problem). 1 


Step 2. 
Step 3. 


P C = P D + i ( p s +P s > As - 

1 1 t Dj S Cl 


The compatibility relation applied along the C+ character- 
istic ACj yields 


Step 4. e 


1 


AS - % 


1 


* y = y - y L where L denotes the lower subsonic boundary. If the lower bound- 
ary is an axis of symmetry, y, =0 and the polynomials take the form 
- .2 . „ .4 , > ... .3 


(y) 


= P Q + P 2 y + P 4 / and e (y) = + ey*. 
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Note, that if the lower subsonic boundary is a wall or axis, the flow de- . 
flection e r is known and no compatibility relation is required. 

4 

Step 5. With the pressure gradient and streamline location 
known, the velocity, temperature and species mass 
fraction may be obtained at (or any mesh point I) 
by using the explicit finite difference formulation 
of Equations (5), (7) and (8). 

Step 6. Knowing P g and at point Equations (AI-1) and AI-2) 

may be inverted to yield (P ) and (e ) . 

Jr y C 

4 4 

Step 7. Since the j polynomial contains three unknowns, and 

two boundary conditions are already known at (0 and e ), 
only one value of 0 at the upper boundary C 2 will be con- 
sistent with both the 0 polynomial and the modified con- 
tinuity equation**. A local iteration to find this value 
of e proceeds as follows: 

Step 7a. A value of e r is assumed. 

4 

Step 7b. The coefficients e Q , e 1 and e 2 at y c , e=e c and 

e=e and at y r , e=e r . The e polynomial is then 

y c x 4 4 

6 (y) = e 0 + 9 1 + S 2 ^ 

ste p 7c - \~--k ( vV'% 2 ' 

(AII-3) 


**(M 2 -1) P s + YPM 2 6 n = A 
where A is described by Equation (AI-3). 
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Step 7d. The compatibility relation along characteristic 
BC 2 yields P c . 


Step 7e 


fp -p ) - p 
AS % H s 


Step 7f . Using the modified continuity equation 


A - (M -1) P, 


9 = - 

n C 

^2 


Step 7g. Since ~ = cose ^ + sine ± ey^ = cose^ 


e„ + sine r e 
\ °2 S C 2 


Step 7h. But using the polynomial j 


e - + 2e 2 y 

l 2 

Step 7i . The values of 0 V and 0y must be identical 

to within a specified tolerance. If they do 
not agree, a new value of is assumed and 
Steps (7b) - (7h) are repeated. Convergence is 
obtained by use of a linear error extrapolation 
procedure . 


Step 8. Flow properties are computed at 


Step 9. P is computed using the normal momentum equation 
C 2 

P n r = ' 6 s r • 

l 2 u 2 

Step 10. P is computed P.. = cose r p„ + sina r P c . 

Jr Jr n r 5 r 
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Step 11* The coefficients of the pressure polynomial between 
and can now be evaluated using the conditions 


at 


y r , P=P r and P =P 

U 1 4 y y ( 


Step 12. 


at y r , P=P C and P =P 

4 4 y y C 2 

yielding P^ y j = p 0 + P^ +■ P^ 2 + P 3 y3 * 

The interior mesh points I between and are evaluated 
as follows. 


Step 12a. A value of yj is assumed. 

Step 12b. Using the polynomial P^ y ^ and e^, Py and 6j 
are obtained. 

Step 12c. Flow properties are obtained using the stream- 
line relations. 

Step 12d. The mass flow at I is evaluated by 

= + Kpq cos0 )i_i + (p 9 cose )|] 

(yj +j - 

1+3 ' 


Step 12e. Since II is a streamline, ipj should equal 
If it does not agree to within a specified 

tolerance, a new value of y^ is assumed and the 
process is iterated until convergence. 
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Step 13. The overall mass flow between is compared to that 
between D^. If the mass flow is not correct the en- 
tire process (Steps (1) - (13)) must then be repeated 
with a new value of P_ for the direct problem. 

V 

L 1 


Generally, convergence was obtained with under five iterations, except when 
the subsonic flow field is close to sonic in a mass averaged sense. This 
situation occurs in the vicinity of the supersonic to subsonic transition 
and again in the vicinity of the subsonic to supersonic transition. In 
these critical regions, the error curve may yield two solutions (supersonic 
and subsonic branch) in which case knowledge of downstream conditions indicates 
which branch is physically correct. For arbitrary initial conditions, the 
flow does not necessarily undergo a smooth super to subsonic transition and 
the matching Mach number tends to be an important factor in effecting this 
transition. The transition from sub to supersonic flow (i.e., closing the 
sonic line) poses an even greater problem. 


In this region, the flow may choke, indicating that a shock is required in 
the initial super to sub transition to allow the mass flow to pass. If 
the subsonic region is bounded by a wall on one side, when choking occurs, 
the wall may be opened to allow passage of the correct mass flow. If the 
subsonic region is embedded in the flow field, then the wall would have to 
be modified upstream of the choked station to allow passage of the proper 
mass flow. In the inverse case, one of the subsonic boundaries would have 
to be a wall and P prescribed at one boundary would yield the wall shape 
at the other boundary. 
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APPENDIX II I 

STEP SIZE AND GRID SPACING CRITERION 

A. Step Size Criterion - Limitations on the marching step ax result 
from characteristic criterion, parabolic stability criterion associated with 
an explicit difference scheme and criterion associated with the solution of 
the finite rate chemistry. The criterion associated with characteristics is 
the Courant-Friedrichs-Lewy condition which states that "the numerical domain 
of dependence of a difference scheme should include the domain of dependence 
of the differential equation or else, convergence is not always possible" (Re- 
ference 9). This criterion is numerically satisfied by intersecting all free 
running characteristics ’emanating from grid points of the supersonic portions 
of the initial profile and obtaining the intersection yielding the minimum 
forward marching step ax ^ as shown in Figure (III-l). 

The stability criterion associated with an explicit parabolic marching scheme 
has been discussed extensively in the literature (Reference 10), based upon 
equations with constant coefficients. The solution of the s-momentum, energy 
and species diffusion equation are all solved by an explicit scheme using the 
grid depicted in Figure ( I I I - 2 ) where the spacing Ay between grid points is 
not necessarily uniform. Linearized stability criterion imposes the condition: 

2y L AX 

7± l (AIII-1) 

pq Pr Re Ay. 

Then the parabolic step size is determined by applying Equation (AIII-1) at each 
grid point I and selecting the minimum ax^. The Ay employed in the equation 
is the average Ay of the interval 1-1 to 1+1; i.e., Ay^ - (Ay^A y 2 ). Then 

P Re 2 

AX vis = FTT Cp I ] . (AIII-2) 

e mm 

It should be noted that the hyperbolic (characteristic) stability criterion is 
determined with the diffusive terms neglected and the parabolic (diffusive) 
stability limit is determined with the- convective terms neglected. It had 
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been thought that the smaller of these two step sizes should be taken as the 
stability limit of the system as noted in Reference (2), but it was also 
pointed out in this reference that use of this criterion still led to in- 
stabilities in some cases. A superposition of the diffusive and convective 
terms may yield more restrictive stability limits than choosing the smaller 
of the step sizes. As pointed out by Cheng in Reference (11), the appropriate 
combined "hyperbol ic-parabol ic“ criterion is 


Ax 


4 


AX 


char 


JL 

+ -L 

ax 


vis 


The finite rate hydrogen-air chemistry employed in this analysis is described 
in References (12) and (13) and uses a fixed time interval At = 4. *10 7 seconds. 
Hence, the chemical step size is 


q T *U 

1 “ i. 

AX . — “ At 

chem r jet 

Since the chemistry is uncoupled based on the procedure developed in Reference 
(14) and described for viscous-characteristics in Reference (2), several chem- 
ical steps may be taken in one overall marching step. The number of chemical 
steps to be allowed in a marching step is left as a user option and is depend- 
ent on the grid spacing in the transverse direction (Ay spacing). Note that 
when mixing is initiated, a very small grid spacing Ay is required, hence the 
overall marching step may be significantly smaller than a chemical step, while 
far downstream, when the Ay is substantially larger, the chemistry criterion 
may dominate the overall step size. 


B. Grid Spacing In The Transverse Direction - The program developed 
can analyze a wide variety of flow situations. The lower boundary may be a 
wall or an axis, while the upper boundary may be a wall or constant pressure 
surface. The initial profile may be entirely nonuniform, in which case all 
grid points must be input initially or may have: (a) uniform jet, uniform 

air (b) uniform jet, nonuniform air (c) nonuniform jet, uniform air. In each 
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of these cases, the jet-air pressure may be balanced or unbalanced. If the 
pressure is unbalanced, the jet underexpansion interaction is calculated. 

For flow with an axis or straight centerbody, only the jet properties, an 
initial Ay based on mixing considerations, and the maximum number of points 
on the jet side have to be input for a uniform jet (cases a or b). As mixing 
progresses, the program adds points at Ay on the jet side, until the maximum 
number of jet points is obtained at which point Ay is doubled, and alternate 
points on the jet side are dropped. This process is continued until the wall 
or axis is reached. Similarly, if the airstream is uniform (case a or c) and 
the upper boundary is such that it sends no waves into the flow field, the 
same procedure applies on the air side. For the unbalanced pressure case, jet 
points are added below the Prandtl -Meyer region and the free stream side of 
the shock serves as the upper boundary if the air is uniform. As the distance 
between the underexpansion shock and the mixing zone increases, mesh points 
are automatically added in this region. 
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APPENDIX IV 
VISCOSITY MODEL 

The turbulent eddy viscosity model referred to as the "Ferri-Kleinstein" model 
is described in References (18) and (19). This model has been incorporated as 
a subroutine in the developed program and assumes an eddy viscosity variation in 
the axial direction only. Incorporation of a different model employing only an 
axial variation simply involves changing the viscosity subroutine while incorpora- 
tion of a model employing both an axial and radial variation can be easily ac- 
complished by minor program modifications. 

Defining x* as the length of the potential core, 

, = Kj X Re ((pu) max - (pU) min ) + K 3 (AIV-1) 

for x < x* 


and 


P = K 2 r h Re ^ <pu) Max - ( P u) min j + K 4 


(AIV-2) 


for x > x* 


where r ? , is the radial distance to the point in the jet where pu = - ((p u ) max + 


^ pu ^min^ * 


For a uniform jet and external stream the term ( P u) max - (pu)^ is replaced by 
the appropriate jet and external stream properties, as depicted in Figure (IV-1). 

For an underexpanded jet, the viscosity is computed on the basis of the region 
below the bounding shock. Due to burning the difference between fp u )j e t and (p u ) ex :'fj 
may be less than that in the combustion zone in which case the max difference is 
obtained as depicted in Figure (IV-2). The coefficients K^, Kg» and K^rand the 
potential core length X* are input items in the program. 
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